library(terra)
ndvi2018 = rast("data/sentinel2_ndvi2018s.tif")
ndvi2023 = rast("data/sentinel2_ndvi2023s.tif")
panel(c(ndvi2018, ndvi2023))
Time for the above code chunk to run: 10.37 seconds
Methods for comparing two layers of spatial continuous raster data
Jakub Nowosad 
This document shows a set of examples related to the article “Comparing spatial patterns in raster data using R”. The examples are focused on comparing two layers of spatial continuous raster data, and are divided into three groups: raster outcome, single value outcome, and multiple value outcome.
The examples in this document use two raster layers representing NDVI values derived from Sentinel-2 images (“Sentinel-2 MSI Level-2A BOA Reflectance” 2022) for early summer of 2018 and 2023 in the area of Tartu, Estonia. The terra package (Hijmans 2024) is used for raster data manipulation.
The calculation time for each example is estimated and displayed at the end of each code chunk.
1 Raster outcome
1.1 Non-spatial context
1.1.1 Non-disjoint areas
1.1.1.1 The difference between values of two rasters for each cell
ndvi_diff = ndvi2023 - ndvi2018
plot(ndvi_diff)
Time for the above code chunk to run: 1.23 seconds
1.2 Spatial context
1.2.1 Non-disjoint areas
1.2.1.1 Correlation coefficient between focal regions of two rasters
ndvi_cor = focalPairs(c(ndvi2023, ndvi2018), w = 5,
fun = "pearson", na.rm = TRUE)
plot(ndvi_cor)
Time for the above code chunk to run: 5.81 seconds
1.2.1.2 The difference between a focal measure of two rasters (1)
This example uses the geodiv package (Smith et al. 2023).
# 'sa': average surface roughness
library(geodiv)
window = matrix(1, nrow = 5, ncol = 5)
ndvi2018_sa_mw = focal_metrics(ndvi2018, window = window,
metric = "sa", progress = FALSE)
ndvi2023_sa_mw = focal_metrics(ndvi2023, window = window,
metric = "sa", progress = FALSE)
ndvi_diff_sa_mw = ndvi2023_sa_mw$sa - ndvi2018_sa_mw$sa
plot(ndvi_diff_sa_mw)
Time for the above code chunk to run: 43.49 seconds
1.2.1.3 The difference between a focal measure of two rasters (2)
This example uses the GLCMTextures package (Ilich 2020).
library(GLCMTextures)
ndvi2018_q = quantize_raster(ndvi2018, n_levels = 16, method = "equal prob")
ndvi2023_q = quantize_raster(ndvi2023, n_levels = 16, method = "equal prob")
ndvi2018_textures = glcm_textures(ndvi2018_q, w = c(5, 5), na.rm = TRUE,
metrics = "glcm_homogeneity",
n_levels = 16, quantization = "none")
ndvi2023_textures = glcm_textures(ndvi2023_q, w = c(5, 5), na.rm = TRUE,
metrics = "glcm_homogeneity",
n_levels = 16, quantization = "none")
ndvi2023_textures_diff = ndvi2023_textures - ndvi2018_textures
plot(ndvi2023_textures_diff)
Time for the above code chunk to run: 211.59 seconds
1.2.1.4 Spatial autocorrelation analysis of the differences (Cliff 1970)
ndvi_diff = ndvi2023 - ndvi2018
ndvi_diff_autocor = autocor(ndvi_diff, method = "moran", global = FALSE)
plot(ndvi_diff_autocor)
Time for the above code chunk to run: 1.02 seconds
1.2.1.5 Structural similarity index (SSIM) between two rasters
This example uses the SSIMmap package (Ha and Long 2023).
library(SSIMmap)
ndvi_ssim = ssim_raster(ndvi2018, ndvi2023, global = FALSE, w = 5)
plot(ndvi_ssim)
Time for the above code chunk to run: 5.5 seconds
1.2.1.6 Comparison of Rao’s quadratic entropy (Rao 1982)
This example uses the rasterdiv package (Rocchini et al. 2021).
library(rasterdiv)
ndvi2018_int = ndvi2018 * 100
ndvi2023_int = ndvi2023 * 100
ndvi2018_rao = paRao(ndvi2018_int, window = 5, progBar = FALSE)
ndvi2023_rao = paRao(ndvi2023_int, window = 5, progBar = FALSE)
ndvi_rao_diff = ndvi2023_rao[[1]][[1]] - ndvi2018_rao[[1]][[1]]
plot(ndvi_rao_diff)
Time for the above code chunk to run: 556.63 seconds
2 Single value outcome
2.1 Non-spatial context
2.1.1 Non-disjoint areas
2.1.1.1 Statistics of the differences between rasters’ values (1)
This example uses the yardstick package (Kuhn, Vaughan, and Hvitfeldt 2024).
library(yardstick)
ndvi_rmse = rmse_vec(values(ndvi2023)[,1], values(ndvi2018)[,1])
ndvi_rmse[1] 0.2191853Time for the above code chunk to run: 0.07 seconds
2.1.1.2 Statistics of the differences between rasters’ values (2)
This example uses the diffeR package (Pontius Jr. and Santacruz 2023).
library(diffeR)
ndvi_mad = MAD(ndvi2023, ndvi2018)[1] "The mean of grid1 is less than the mean of grid2"ndvi_mad Resolution Quantity Strata Element Total
1 10 0.1783152 0 0.005990721 0.184306Time for the above code chunk to run: 0.36 seconds
2.1.2 Disjoint areas
2.1.2.1 Disimilarity between the distributions of two rasters’ values
This example uses the philentropy (HG 2018) package.
library(philentropy)
softmax = function(x) {
exp_x = exp(x - max(x))
return(exp_x / sum(exp_x))
}
ndvi2018_vals = na.omit(values(ndvi2018)[,1])
ndvi2023_vals = na.omit(values(ndvi2023)[,1])
ndvi2018_vals_prob = softmax(ndvi2018_vals)
ndvi2023_vals_prob = softmax(ndvi2023_vals)
ndvi2023_vals_prob_interp = approx(seq_along(ndvi2023_vals_prob),
ndvi2023_vals_prob,
xout = seq_along(ndvi2018_vals_prob))$y
ndvi_mat = rbind(ndvi2023_vals_prob_interp, ndvi2018_vals_prob)
philentropy::distance(ndvi_mat, method = "kullback-leibler")kullback-leibler
0.03759518 Time for the above code chunk to run: 0.11 seconds
2.2 Spatial context
2.2.1 Non-disjoint areas
2.2.2 Disjoint areas
2.2.2.1 The difference between an average of a focal measure of two rasters (1)
This example uses the geodiv package (Smith et al. 2023).
library(geodiv)
ndvi2018_sa = sa(ndvi2018)
ndvi2023_sa = sa(ndvi2023)
abs(ndvi2023_sa - ndvi2018_sa)[1] 0.06886273Time for the above code chunk to run: 0.05 seconds
2.2.2.2 The difference between an average of a focal measure of two rasters (1)
This example uses the GLCMTextures package (Ilich 2020).
library(GLCMTextures)
ndvi2018_q = quantize_raster(ndvi2018, n_levels = 16, method = "equal prob")
ndvi2023_q = quantize_raster(ndvi2023, n_levels = 16, method = "equal prob")
ndvi2018_q_mat = as.matrix(ndvi2018_q, wide = TRUE)
ndvi2023_q_mat = as.matrix(ndvi2023_q, wide = TRUE)
ndvi2018_horizontal_glcm = make_glcm(ndvi2018_q_mat, n_levels = 16,
shift = c(1,0), na.rm = TRUE)
ndvi2023_horizontal_glcm = make_glcm(ndvi2023_q_mat, n_levels = 16,
shift = c(1,0), na.rm = TRUE)
ndvi2018_hom = glcm_metrics(ndvi2018_horizontal_glcm, "glcm_homogeneity")
ndvi2023_hom = glcm_metrics(ndvi2023_horizontal_glcm, "glcm_homogeneity")
abs(ndvi2023_hom - ndvi2018_hom)glcm_homogeneity
0.01880257 Time for the above code chunk to run: 0.5 seconds
2.2.2.3 Average of Structural Similarity Index (Wang et al. 2004; Robertson et al. 2014)
This example uses the SSIMmap package (Ha and Long 2023).
library(SSIMmap)
ssim_raster(ndvi2018, ndvi2023, global = TRUE)SSIM: 0.63845 SIM: 0.90432 SIV: 0.90778 SIP: 0.75671 Time for the above code chunk to run: 0.85 seconds
2.2.2.4 The similarity index (CQ) based on the co-dispersion coefficient
This example uses the SpatialPack package (Vallejos, Osorio, and Bevilacqua 2020).
library(SpatialPack)
ndvi2018_mat = as.matrix(ndvi2018, wide = TRUE)
ndvi2023_mat = as.matrix(ndvi2023, wide = TRUE)
ndvi_CQ = CQ(ndvi2018_mat, ndvi2023_mat)
ndvi_CQ$CQTime for the above code chunk to run: 0 seconds
2.2.2.5 Comparison of the values of the Boltzmann entropy of a landscape gradient (Gao and Li 2019)
This example uses the bespatial package (Nowosad 2024).
library(bespatial)
ndvi2018_bes = bes_g_gao(ndvi2018, method = "hierarchy", relative = TRUE)
ndvi2023_bes = bes_g_gao(ndvi2023, method = "hierarchy", relative = TRUE)
abs(ndvi2023_bes$value - ndvi2018_bes$value)[1] 15402.11Time for the above code chunk to run: 0.18 seconds
2.2.2.6 Comparison of deep learning-based feature maps using a dissimilarity measure (Malik and Robertson 2021)
This example uses the keras3 (Kalinowski, Allaire, and Chollet 2024) and philentropy (HG 2018) packages.
library(keras3)
library(philentropy)
# keras3::install_keras(backend = "tensorflow")
normalize_raster = function(r) {
min_val = terra::global(r, "min", na.rm = TRUE)[[1]]
max_val = terra::global(r, "max", na.rm = TRUE)[[1]]
r = terra::app(r, fun = function(x) (x - min_val) / (max_val - min_val))
return(r)
}
ndvi_2018n = normalize_raster(ndvi2018)
ndvi_2023n = normalize_raster(ndvi2023)
ndvi2018_mat = as.matrix(ndvi_2018n, wide = TRUE)
ndvi2023_mat = as.matrix(ndvi_2023n, wide = TRUE)
ndvi2018_mat = array(rep(ndvi2018_mat, 3),
dim = c(nrow(ndvi2018_mat), ncol(ndvi2018_mat), 3))
ndvi2023_mat = array(rep(ndvi2023_mat, 3),
dim = c(nrow(ndvi2023_mat), ncol(ndvi2023_mat), 3))
model = keras3::application_vgg16(weights = "imagenet", include_top = FALSE,
input_shape = c(nrow(ndvi2018_mat), ncol(ndvi2018_mat), 3))
ndvi2018_mat = keras3::array_reshape(ndvi2018_mat, c(1, dim(ndvi2018_mat)))
ndvi2023_mat = keras3::array_reshape(ndvi2023_mat, c(1, dim(ndvi2023_mat)))
features2018 = predict(model, ndvi2018_mat)1/1 - 3s - 3s/stepfeatures2023 = predict(model, ndvi2023_mat)1/1 - 2s - 2s/stepfeature_map_2018_1 = as.vector(features2018[1,,,1])
feature_map_2023_1 = as.vector(features2023[1,,,1])
distance(rbind(feature_map_2018_1, feature_map_2023_1))euclidean
3.446775 Time for the above code chunk to run: 15.43 seconds
3 Multiple values outcome
3.0.0.1 The distribution of the difference between values of two rasters
ndvi_diff = ndvi2023 - ndvi2018
hist(ndvi_diff)
Time for the above code chunk to run: 0.15 seconds
3.0.0.2 Statistics of the differences between rasters’ values calculated at many scales
This example uses the waywiser package (Mahoney 2023).
library(waywiser)
cell_sizes = c(200, 400, 600)
ndvi_multi_scale = ww_multi_scale(truth = ndvi2018, estimate = ndvi2023,
metrics = list(yardstick::rmse),
cellsize = cell_sizes,
progress = FALSE)
ndvi_multi_scale# A tibble: 3 × 6
.metric .estimator .estimate .grid_args .grid .notes
<chr> <chr> <dbl> <list> <list> <list>
1 rmse standard 0.193 <tibble [1 × 1]> <sf [625 × 5]> <tibble [0 × 2]>
2 rmse standard 0.191 <tibble [1 × 1]> <sf [169 × 5]> <tibble [0 × 2]>
3 rmse standard 0.188 <tibble [1 × 1]> <sf [81 × 5]> <tibble [0 × 2]>Time for the above code chunk to run: 3.47 seconds
References
Cliff, Andrew D. 1970. “Computing the Spatial Correspondence Between Geographical Patterns.” Transactions of the Institute of British Geographers, no. 50 (July): 143. https://doi.org/10.2307/621351.
Gao, Peichao, and Zhilin Li. 2019. “Aggregation-Based Method for Computing Absolute Boltzmann Entropy of Landscape Gradient with Full Thermodynamic Consistency.” Landscape Ecology 34 (8): 1837–47. https://doi.org/gkb47k.
Ha, Hui Jeong (Hailyee), and Jed Long. 2023. SSIMmap: The Structural Similarity Index Measure for Maps. Manual.
HG, Drost. 2018. “Philentropy: Information Theory and Distance Quantification with R.” Journal of Open Source Software 3 (26): 765.
Hijmans, Robert J. 2024. Terra: Spatial Data Analysis. Manual.
Ilich, Alexander R. 2020. “GLCMTextures.” https://doi.org/10.5281/zenodo.4310186.
Kalinowski, Tomasz, JJ Allaire, and François Chollet. 2024. Keras3: R Interface to ’Keras’. Manual.
Kuhn, Max, Davis Vaughan, and Emil Hvitfeldt. 2024. Yardstick: Tidy Characterizations of Model Performance. Manual.
Mahoney, Michael J. 2023. “Waywiser: Ergonomic Methods for Assessing Spatial Models.” https://doi.org/10.48550/arXiv.2303.11312.
Malik, Karim, and Colin Robertson. 2021. “Landscape Similarity Analysis Using Texture Encoded Deep-Learning Features on Unclassified Remote Sensing Imagery.” Remote Sensing 13 (3): 492. https://doi.org/10.3390/rs13030492.
Nowosad, Jakub. 2024. Bespatial: Boltzmann Entropy for Spatial Data. Manual.
Pontius Jr., Robert Gilmore, and Ali Santacruz. 2023. diffeR: Metrics of Difference for Comparing Pairs of Maps or Pairs of Variables. Manual.
Rao, C.Radhakrishna. 1982. “Diversity and Dissimilarity Coefficients: A Unified Approach.” Theoretical Population Biology 21 (1): 24–43. https://doi.org/10.1016/0040-5809(82)90004-1.
Robertson, Colin, Jed A. Long, Farouk S. Nathoo, Trisalyn A. Nelson, and Cameron C. F. Plouffe. 2014. “Assessing Quality of Spatial Models Using the Structural Similarity Index and Posterior Predictive Checks.” Geographical Analysis 46 (1): 53–74. https://doi.org/10.1111/gean.12028.
Rocchini, Duccio, Elisa Thouverai, Matteo Marcantonio, Martina Iannacito, Daniele Da Re, Michele Torresani, Giovanni Bacaro, et al. 2021. “Rasterdiv - An Information Theory Tailored R Package for Measuring Ecosystem Heterogeneity from Space: To the Origin and Back.” Methods in Ecology and Evolution 12 (6): 2195. https://doi.org/10.1111/2041-210X.13583.
“Sentinel-2 MSI Level-2A BOA Reflectance.” 2022. European Space Agency. https://doi.org/10.5270/s2_-znk9xsj.
Smith, Annie C., Phoebe Zarnetske, Kyla Dahlin, Adam Wilson, and Andrew Latimer. 2023. Geodiv: Methods for Calculating Gradient Surface Metrics. Manual.
Vallejos, Ronny, Felipe Osorio, and Moreno Bevilacqua. 2020. Spatial Relationships Between Two Georeferenced Variables: With Applications in R. Springer Nature.
Wang, Z., A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli. 2004. “Image Quality Assessment: From Error Visibility to Structural Similarity.” IEEE Transactions on Image Processing 13 (4): 600–612. https://doi.org/10.1109/TIP.2003.819861.